WHEN CORRELATIONS MATTER - RESPONSE OF DYNAMICAL NETWORKS TO 

SMALL PERTURBATIONS 



Thimo Rohlf 1,2 , Natali Gulbahce 3,4 and ChristofTeuscher 5 

x Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA 
2 Max-Planck-Institute for Mathematics in the Sciences, Inselstrasse 22, D-04103 Leipzig, Germany 
3 Center for Complex Network Research, Northeastern University, Boston, MA 021 15, USA 
4 Center for Cancer Systems Biology, Dana Farber Cancer Institute, Boston, MA, 02215, USA 
5 Los Alamos National Laboratory, CCS-3, MS B256, Los Alamos, NM 87545, USA 

rohlf @ santafe.edu 



ABSTRACT 

We systematically study and compare damage spreading 
for random Boolean and threshold networks under small 
external perturbations (damage), a problem which is rele- 
vant to many biological networks. We identify a new char- 
acteristic connectivity K s , at which the average number of 
damaged nodes after a large number of dynamical updates 
is independent of the total number of nodes N, We esti- 
mate the critical connectivity for finite N and show that it 
systematically deviates from the annealed approximation. 
Extending the approach followed in a previous study ITTIl . 
we present new results indicating that internal dynamical 
correlations tend to increase not only the probability for 
small, but also for very large damage events, leading to a 
broad, fat-tailed distribution of damage sizes. These find- 
ings indicate that the descriptive and predictive value of 
averaged order parameters for finite size networks - even 
for biologically highly relevant sizes up to several thou- 
sand nodes - is limited. 

1. INTRODUCTION 

Random Boolean networks (RBN) were originally intro- 
duced as simplified models of gene regulation [5, 6|. In 
the limit of large system sizes, they exhibit a dynami- 
cal order-disorder transition at a critical wiring density 
K c (4); similar observations were made for sparsely con- 
nected random threshold (neural) networks (RTN) ll7l [T0ll . 
For a finite system size N, the dynamics of both systems 
converge to periodic attractors after a finite number of up- 
dates. At K c , the phase space structure in terms of attrac- 
tor periods [ 1 ], the number of different attractors [13] and 
the distribution of basins of attraction f2) is complex. Fur- 
thermore, critical networks exhibit many properties rem- 
iniscent of biological networks, leading to the idea K c 
might be an "attractor of evolution" 0. 

To ensure proper function, regulatory networks in liv- 
ing cells have to be robust (insensitive) against external 
perturbations. In terms of RBN/RTN dynamics, perturba- 
tions can disrupt the generic dynamical state (fixed point 
or periodic attractor) of the network, and hence are re- 
ferred to as "damage"; this type of study has been applied, 



for example, to the perturbation of gene expression pat- 
terns in a cell due to mutations . 

Mean-field techniques as, for example, the annealed 
approximation (AA) introduced by Derrida and Pomeau 
H, allow for an analytical treatment of damage spreading 
and exact determination of the critical connectivity K c un- 
der various constraints [14]. It has been shown that local 
rewiring rules coupled to mean-field-like order parame- 
ters of the dynamics can drive both RBN and RTN to self- 
organized criticality J3] 0. 

Studies of RBN/RTN dynamics based on the AA usu- 
ally implicitly assume that, at least for large N, principal 
properties of damage spreading should not depend on the 
initial perturbation size. For example, the determination 
of K c using a one-bit initial perturbation (sparse percola- 
tion limit), or an initial perturbation size increasing with 
N should yield the same value for large N, since it is as- 
sumed that correlations can be neglected in this limit by 
averaging over a large number of different random net- 
work realizations. In this paper, we extend results of a 
previous study ifTTI and present the following findings that 
are, at least in part, in clear contradiction to these assump- 
tions: 

• In section 3. 1, we identify a new characteristic point 
K s < K c , where the expectation value of the num- 
ber of damaged nodes after large number of dynam- 
ical updates is independent of N. 

• By the definition of marginal damage spreading, we 
estimate the critical connectivity K C (N) for finite 
N, and present evidence that, even in the large N 
limit, for small initial perturbations K c systemati- 
cally deviates from the predictions of the AA (sec- 
tion 3.2). 

• In section 3.3, we present new results proving that, 
slightly below K c , starting from random initial con- 
ditions, the AA holds only for small times t, indi- 
cating that after passing transient dynamics inher- 
ent correlations considerably affect damage propa- 
gation. 
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Figure 1 . Average Hamming distance (damage) d after 200 sys- 
tem updates, averaged over 10000 randomly generated networks 
for each value of K, with 100 different random initial condi- 
tions and one-bit perturbed neighbor configurations for each net- 
work. For both RBN and RTN, all curves for different TV approx- 
imately intersect in a characteristic point K s . 

• Last, we show that vanishing, as well as large dam- 
age events are overrepresented in damage size statis- 
tics, leading to highly skewed distributions, which 
are poorly characterized by averages (section 3.4). 

2. DYNAMICS 

2.1. Random Boolean Networks 

A RBN is a discrete dynamical system composed of N 
automata. Each automaton is a Boolean variable with two 
possible states: {0, 1}, and the dynamics is such that 

F: {0,1}^ ^{0,1}*, (1) 

where F = (/i, /j, /at), and each fa is represented 
by a look-up table of K inputs randomly chosen from the 
set of N automata. Initially, Ki neighbors and a look-table 
are assigned to each automaton at random. 

An automaton state a\ € {0, 1} is updated using its 
corresponding Boolean function: 

= fM x ,A^-A K )- (2) 

We randomly initialize the states of the automata (initial 
condition of the RBN). The automata are updated syn- 
chronously using their corresponding Boolean functions. 

2.2. Random Threshold Networks 

An RTN consists of A" randomly interconnected binary 
sites (spins) with states cr, = ±1. For each site i, its state 
at time t + 1 is a function of the inputs it receives from 
other spins at time t: 

a,{t + 1) = sgn c V a o (*) + h - j < 3 > 

The N network sites are updated synchronously. In the 
following discussion the threshold parameter h is set to 
zero. The interaction weights Cy take discrete values Cy = 
+ 1 or —1 with equal probability. If i does not receive sig- 
nals from j, one has = 0. 
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Figure 2. The critical connectivity K a c parae (N) in the SP limit 
as a function of N. Curves are power-law fits according to Eq. 
©, straight dashed lines mark K^ nnealed and K s for RBN and 
RTN, respectively. 

3. RESULTS 

3.1. Scaling 

We first study the expectation value d of damage, quan- 
tified by the Hamming distance of two different system 
configurations, after a large number T of system updates. 
Fig. Q]shows d as a function of the average connectivity K 
for different network sizes N by using a random ensem- 
ble for statistics. For both RBN and RTN, the observed 
functional behavior strongly suggests that the curves ap- 
proximately intersect at a common point (K s , d s ), where 
the observed Hamming distance for large t is independent 
of the system size N. 

We verified this finding quantitatively by using finite- 
size-scaling methods ifTTI . In particular, one can show that 
d as a function of N and K obeys the following scaling 
ansatz: 

d(K, N) = a{K) ■ N l{k) + d (K), -1 < 7 < 1. (4) 

It is straight-forward to show that 7^—1 for small 
K — > 0, and that 7 — * 1 for densely connected networks 
above the percolation transition (K > K c ). Evidently, 
this implies that at some characteristic connectivity K s , 
there has to be a transition from negative to positive 7 val- 
ues, with ^{Kg) w 0. It is a very interesting question 
whether K s coincides with K c , or if it is different from 
K c for large N. For a precise numerical determination of 
K s , one can make use of the fact that d exhibits an expo- 
nential dependence near K c : 

d{K, N) m ci (N) exp [c 2 (N) N a K] (5) 

with a w 0.42. High-accuracy fits of this dependence 
(with ci and c 2 as adjustable parameters) in the interval 
1.6 < K < 2.1 yield 

{K* BN ,df BN ) = (1.875 ±0.05, 0.62 ±0.05) (6) 

for RBN and, correspondingly, 

{Kf TN ,df TN ) = (1.729 ± 0.045, 0.51 ± 0.04) (7) 

for RTN. We verified these findings up to N = 16384, 
waiting T = 5000 updates for the dynamics to relax; 
for even larger N, simulations become intractable due to 



exponentially increasing relaxation times. Evidently, we 
tend to miss large damage events since they need the most 
time to develop. Facing this unavoidable biased under- 
sampling of large avalanches, one can argue that the true 
values of K s are probably even lower than our measured 
values. From this evidence, and also from more refined 
scaling arguments [11], we conclude that K s is distinct 
from K c in the limit of large N. 

3.2. Deviations of K c from the annealed approxima- 
tion 

Interestingly, K s is close to, but distinct from the criti- 
cal connectivities K* BN = 2 and K* TN = 1.845, as 
predicted by the AA. Since in this study we consider the 
limit of very weak initial perturbations which is usually 
not covered in theoretical studies of RBN/RTN dynam- 
ics, we now have to consider the possibility that K c itself 
may deviate from the prediction of the AA. An intuitive 
definition of criticality for finite N can be formulated in 
terms of marginal damage spreading. If at time t one bit 
is flipped, one requires at time t+1 lfT4l[T0ll 

d(t + l) = {p s )(K c )K c =l, (8) 

where (p s ) (K) is the average damage propagation proba- 
bility. Fig. 2 shows K^ parse (N), using the values ci(N) 
and C2 ( N) obtained from numerical fits of Eq. for both 
RBN and RTN. We find that both systems, in a very good 
approximation, obey the scaling relationship 

K* parse (N) b ■ N~ s + (9) 

with b = 3.27 ± 0.79, 5 = 0.85 ± 0.07 and = 
1.9082 ± 0.008 for RBN and b = 3.853 ± 0.76, 5 = 
0.736 ±0.05 and Kf = 1.7595 ±0.008 for RTN. Hence, 
in the limit N — > oo, we can extrapolate 

R oo,RBN = 1 90g2 ± o QQg ( 10 ) 

for RBN, and for RTN 

K™> RTN = 1.7595 ± 0.008. (11) 

Thus, for both RBN and RTN in the sparse percolation 
limit, we make the surprising observation that K sparse 
systematically deviates from x^ nnealed . While we find 
K*P arse (N) > K* nnealed for small N < 128, for larger 
N we observe a monotonic decay that approaches an asymp- 
totic value considerably below X^ nneal , suggesting that 
the observed deviations from the AA also hold in the large 
N limit. In the following two subsections, we will extend 
this analysis and discuss possible causes for these devia- 
tions. 

3.3. Time dependence of d 

Since we found systematic deviations from the AA for 
large t, it is interesting to ask whether the AA still holds 
for small f, starting from random initial states. In particu- 
lar, one can derive the following recursive map for damage 
propagation at t > [ 10 1: 

d(t) = N ■ ( Ps )(K) ■ (l - e -K<t-x)/N^ ) (12) 




Figure 3. Time-dependence of (average) damage propagation in 
RTN of size N = 4096 just below K c \ damage d at time t was 
averaged over 10 network realizations and 100 different ini- 
tial conditions (and the corresponding neighbor states with one 
bit perturbed at random) at t = for each data point . Lined 
curves are the corresponding solutions of the AA (Eq. d!2t). For 
t > 20, pronounced deviations of simulation results from the 
AA are found, in particular for K — 1.8. Arrows indicate re- 
sults with "corrected" statistics, i.e. without "pseudo-damage" 
due to attractor phase lags. 

where ( Ps )(K) is the average probability that a link propa- 
gates damage. Let us now test this relationship in the inter- 
esting range K s < K < K^ nnealed for ensembles of ran- 
domly generated networks (RTN with Poissonian degree- 
distribution), with one-bit perturbations of randomly cho- 
sen initial conditions. Figure 3 shows that, for small t, 
the dependence for d(t) found in numerical simulations 
obeys this prediction very well. However, after an initial 
decrease of d(t), an increase above the initial damage size 
(i.e. supercritical behavior) is found, in clear contradiction 
to the AA. This indicates that, after the system has passed 
transient dynamics, inherent dynamical correlations con- 
siderably modify damage propagation (fractal structure of 
attraction basins J2j). One can also show that "pseudo- 
damage" events , i.e. cases where networks run on the 
same attractor, but with a phase lag captured in a non-zero 
Hamming distance, do not substantially contribute to this 
effect (arrows in Fig. 3). This proves that our results are 
very robust against changes in the way statistics is taken. 

3.4. Distribution of damage sizes 

Let us now go beyond averaged (mean-field) quantities 
and investigate detailed statistics of damage sizes. For this 
purpose, for different K and ensembles of Z e random 
network realizations were created; for each network real- 
ization, Zi random initial conditions a (plus a neighbor 
state with one bit perturbed at random) were tested, and 
statistics of damage sizes was taken after 1000 dynami- 
cal updates. Notice that we do not average damage sizes 
for a given network realization, since this would again 
represent a kind of mean-field approximation. Figure 4 
shows that the resulting statistical distributions near K c 
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Figure 4. Statistical distribution p(d) of damage sizes for three 
different system sizes: N = 64 (+), N = 256 (x) and N = 
1024 (*). Lined curves are solutions of Eq. | |13I >. 



are highly skewed, with more than 90% events of van- 
ishing damage size, and a fiat tail of large damage events 
which becomes more and more pronounced for increasing 
N. Similar problems have been studied by Samuelsson 
and Socolar [ 12 1 for the number of undamaged nodes u in 
the limit of exhaustive percolation. From symmetry con- 
siderations, it follows that the probability distribution p(d) 
of the number d of damaged nodes in the limit of sparse 
percolation obeys a similar dependence as u in the case of 
exhaustive percolation, and hence 



p(d) « a(N) 



exp| 



\{d-N~ 2 l*y 



Vd • TV- 2 / 3 



(13) 



where a(N) is a free parameter. One finds that the results 
of numerical simulations agree very well with this esti- 
mate even for considerably small N (Fig. 4). From the 
shape of these distributions, one recognizes that vanish- 
ing, as well as large damage events are much more proba- 
ble than expected from mean-field considerations. In part, 
this explains the deviations from the annealed approxima- 
tion found for d near criticality (Fig. 3), and it also ques- 
tions in how far averaged quantities deliver an informative 
description of RBN/RTN dynamics for finite size N. 

4. DISCUSSION 

We showed that, for very weak (one-bit) perturbations of 
the initial states of RBN and RTN dynamics, the result- 
ing damage at later times exhibits a non-trivial scaling 
with network size N, and, near the critical order-disorder 
transition - the so-called the 'edge of chaos' - consider- 
able deviations from the annealed approximation. These 
deviations have escaped earlier studies, since usually the 
rescaled damage d/N (or the overlap 1 — d/N, respec- 
tively) was studied, and the thermodynamic limit of large 
N was considered. Our study indicates that there is a 
strong need for more refined studies of damage propaga- 
tion in RBN/RTN, that explicitly take into account dy- 
namical correlations and the fractal structure of attrac- 
tion basins [2J. One may expect that the situation is even 
more complex for networks with more realistic topolo- 
gies. Even for simple random graphs, as applied in this 



study, damage size distributions are highly skewed, ques- 
tioning the descriptive and predictive value of simple, av- 
eraged order parameters for this class of complex systems. 
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